#
# data <- readRDS('data_sim_p_1000.rds')
# for(x in names(data)) assign(x, data[[x]], envir = globalenv())
# U_full <- data$U
# rm(data)
#
# S.data.time <- S.data$obs
# S.data.time.log.sum <- sum(log(S.data.time))
# Exemple d'optimalité entre les 3 fct MH, MH_high_dim et MH_high_dim_para sur le modèle non lineaire mixte
m <- function(t, phi1, phi2, phi3) (phi1 )/(1+exp((phi2-t)/phi3))
#=======================================#
p <- 1000
parameter <- list(sigma2 = .05^2,
mu = c(0.9,90,5),
omega2 = c(0.005, 40, 1),
bara = 90,
barb = 30,
baralpha = 0.5,
beta = rep(0,p))
parameter$beta[1:4] <- c(-.8, -.2 , .3 , .9)
#=======================================#
t <- seq(60,120, length.out = 10) #time values
set.seed(123)
G <- 40 ; ng = 4
link = m
dt <- create_JLS_HD_data(G, ng, t, m, link, parameter)
var.true <- dt$var.true
a <- var.true$a ; var.true$a <- NULL #a fixé (et retiré des variables latentes)
S.data <- dt$survival
U_full <- dt$U
Y <- do.call(get_obs, var.true) + rnorm(n, 0, sqrt(parameter$sigma2))
S.data.time <- S.data$obs
S.data.time.log.sum <- sum(log(S.data.time))
model <- SAEM_model(
function(sigma2, ...) -n/(2*sigma2),
function(phi1, phi2, phi3, ...) mean((Y - get_obs(phi1, phi2, phi3) )^2 ), 'sigma2',
# === Variable Latente === #
latent_vars = list(
# === Non linear model === #
latent_variable('phi', dim = G, size = 3, prior = list(mean = 'mu', variance = 'omega2'),
add_on = c('zeta(phi1 = phi1, phi2 = phi2, phi3 = phi3, ...)' )),
# === S.data model === #
latent_variable('b', prior = list(mean = 'barb', variance.hyper = 'sigma2_b'),
add_on = c('zeta(b = b, ...) +',
'sum(h$eval(b = b, ..., i = c(1,2)))' )),
latent_variable('alpha', prior = list(mean = 'baralpha', variance.hyper = 'sigma2_alpha'),
add_on = c('zeta(alpha = alpha, ...) +',
'alpha*h$eval(alpha = alpha,..., i = 3)'))
),
# === Paramètre de regression === #
regression.parameter = list(
regression_parameter('beta', 1, function(...) SPGD(5, theta0 = beta,
step = 0.05, lambda = 1/sqrt(N),
normalized.grad = T,
zeta.der.B, N, zeta.B,
Z$alpha, Z$phi1, Z$phi2, Z$phi3,Z$b) )
)
)
# --- Initialisation des paramètres --- #
parameter0 <- parameter %>% sapply(function(x) x* runif(1, 1.1,1.4))
#===============================================#
load.SAEM(model)
U <- U_full
oracle <- maximisation(1, do.call(S$eval, var.true), parameter, var.true)
#==============================================================================#
init.options <- list(x0 = list(phi = c(1,80,4), b = parameter0$barb, alpha = parameter0$baralpha),
sd = list(phi = c(.05, 1.5, .5), b = 1, alpha = .1) )
SAEM.options <- list(niter = 300, sim.iter = 5, burnin = 190,
adptative.sd = 0.6)
saem
p <- 4
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 48min 37sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9033
|
89.9108
|
5.0104
|
0.0040
|
35.2265
|
0.5900
|
29.6106
|
0.7293
|
|
Rrmse
|
0.0022
|
0.0001
|
0.0005
|
0.0005
|
0.0105
|
0.0192
|
0.1507
|
0.0130
|
0.4585
|
p <- 20
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 48min 51sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9033
|
89.9168
|
5.0148
|
0.0040
|
35.1852
|
0.5881
|
28.5855
|
0.8915
|
|
Rrmse
|
0.0025
|
0.0002
|
0.0005
|
0.0014
|
0.0207
|
0.0204
|
0.1535
|
0.0471
|
0.7830
|
p <- 100
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 54min 04sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9035
|
89.9555
|
5.0277
|
0.0039
|
35.2423
|
0.6457
|
44.0060
|
2.1571
|
|
Rrmse
|
0.0006
|
0.0003
|
0.0000
|
0.0040
|
0.0095
|
0.0188
|
0.0706
|
0.4669
|
3.3142
|
p <- 500
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 55min 16sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9034
|
89.9250
|
5.0230
|
0.0039
|
35.5875
|
0.5844
|
35.969
|
1.3529
|
|
Rrmse
|
0.0017
|
0.0003
|
0.0004
|
0.0030
|
0.0007
|
0.0092
|
0.1588
|
0.199
|
1.7059
|
p <- 1000
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 54min 29sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
0.5000
|
|
Estimated value
|
0.0025
|
0.9039
|
89.9224
|
5.0250
|
0.0039
|
35.3015
|
0.5846
|
30.6821
|
0.6981
|
|
Rrmse
|
0.0018
|
0.0008
|
0.0004
|
0.0034
|
0.0026
|
0.0172
|
0.1585
|
0.0227
|
0.3963
|